/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | foam-extend: Open Source CFD
   \\    /   O peration     | Version:     4.1
    \\  /    A nd           | Web:         http://www.foam-extend.org
     \\/     M anipulation  | For copyright notice see file Copyright
-------------------------------------------------------------------------------
License
	This file is part of foam-extend.

	foam-extend is free software: you can redistribute it and/or modify it
	under the terms of the GNU General Public License as published by the
	Free Software Foundation, either version 3 of the License, or (at your
	option) any later version.

	foam-extend is distributed in the hope that it will be useful, but
	WITHOUT ANY WARRANTY; without even the implied warranty of
	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
	General Public License for more details.

	You should have received a copy of the GNU General Public License
	along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.

Description

\*---------------------------------------------------------------------------*/

#include "surfaceSets.H"
#include "polyMesh.H"
#include "triSurface.H"
#include "triSurfaceSearch.H"
#include "pointSet.H"
#include "cellSet.H"
#include "surfaceToCell.H"
#include "cellToPoint.H"
#include "cellToCell.H"
#include "pointToCell.H"
#include "meshSearch.H"
#include "cellClassification.H"


// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //


// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //

//Foam::scalar Foam::surfaceSets::minEdgeLen
//(
//    const primitiveMesh& mesh,
//    const label pointI
//)
//{
//    const edgeList& edges = mesh.edges();
//
//    const pointField& points = mesh.points();
//
//    const labelList& pEdges = mesh.pointEdges()[pointI];
//
//    scalar minLen = GREAT;
//
//    forAll(pEdges, i)
//    {
//        minLen = min(minLen, edges[pEdges[i]].mag(points));
//    }
//    return minLen;
//}
//
//
//// Returns true if cell uses at least one selected point
//bool Foam::surfaceSets::usesPoint
//(
//    const primitiveMesh& mesh,
//    const boolList& selectedPoint,
//    const label cellI
//)
//{
//    const labelList& cFaces = mesh.cells()[cellI];
//
//    forAll(cFaces, cFaceI)
//    {
//        label faceI = cFaces[cFaceI];
//
//        const face& f = mesh.faces()[faceI];
//
//        forAll(f, fp)
//        {
//            if (selectedPoint[f[fp]])
//            {
//                return true;
//            }
//        }
//    }
//    return false;
//}



//// Remove cells in allCells which are connected to other cells in allCells
//// by outside vertices only. Since these outside vertices will be moved onto
//// a surface they might result in flat cells.
//Foam::label Foam::surfaceSets::removeHangingCells
//(
//    const primitiveMesh& mesh,
//    const triSurfaceSearch& querySurf,
//    labelHashSet& internalCells
//)
//{
//    const pointField& points = mesh.points();
//    const cellList& cells = mesh.cells();
//    const faceList& faces = mesh.faces();
//
//    // Determine cells that have all points on the boundary.
//    labelHashSet flatCandidates(getHangingCells(mesh, internalCells));
//
//    // All boundary points will become visible after subsetting and will be
//    // snapped
//    // to surface. Calculate new volume for cells using only these points and
//    // check if it does not become too small.
//
//    // Get points used by flatCandidates.
//    labelHashSet outsidePoints(flatCandidates.size());
//
//    // Snap outside points to surface
//    pointField newPoints(points);
//
//    for
//    (
//        labelHashSet::const_iterator iter = flatCandidates.begin();
//        iter != flatCandidates.end();
//        ++iter
//    )
//    {
//        const cell& cFaces = cells[iter.key()];
//
//        forAll(cFaces, cFaceI)
//        {
//            const face& f = faces[cFaces[cFaceI]];
//
//            forAll(f, fp)
//            {
//                label pointI = f[fp];
//
//                if (outsidePoints.insert(pointI))
//                {
//                    // Calculate new position for this outside point
//                    tmp<pointField> tnearest =
//                        querySurf.calcNearest(pointField(1, points[pointI]));
//                    newPoints[pointI] = tnearest()[0];
//                }
//            }
//        }
//    }
//
//
//    // Calculate new volume for mixed cells
//    label nRemoved = 0;
//    for
//    (
//        labelHashSet::const_iterator iter = flatCandidates.begin();
//        iter != flatCandidates.end();
//        ++iter
//    )
//    {
//        label cellI = iter.key();
//
//        const cell& cll = cells[cellI];
//
//        scalar newVol = cll.mag(newPoints, faces);
//        scalar oldVol = mesh.cellVolumes()[cellI];
//
//        if (newVol < 0.1 * oldVol)
//        {
//            internalCells.erase(cellI);
//            nRemoved++;
//        }
//    }
//
//    return nRemoved;
//}


//// Select all points out of pointSet where the distance to the surface is less
//// than a factor times a local length scale (minimum length of connected edges)
//void Foam::surfaceSets::getNearPoints
//(
//    const primitiveMesh& mesh,
//    const triSurface&,
//    const triSurfaceSearch& querySurf,
//    const scalar edgeFactor,
//    const pointSet& candidateSet,
//    pointSet& nearPointSet
//)
//{
//    if (edgeFactor <= 0)
//    {
//        return;
//    }
//
//    labelList candidates(candidateSet.toc());
//
//    pointField candidatePoints(candidates.size());
//    forAll(candidates, i)
//    {
//        candidatePoints[i] = mesh.points()[candidates[i]];
//    }
//
//    tmp<pointField> tnearest = querySurf.calcNearest(candidatePoints);
//    const pointField& nearest = tnearest();
//
//    const pointField& points = mesh.points();
//
//    forAll(candidates, i)
//    {
//        label pointI = candidates[i];
//
//        scalar minLen = minEdgeLen(mesh, pointI);
//
//        scalar dist = mag(nearest[i] - points[pointI]);
//
//        if (dist < edgeFactor * minLen)
//        {
//            nearPointSet.insert(pointI);
//        }
//    }
//}


// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //


// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

void Foam::surfaceSets::getSurfaceSets
(
	const polyMesh& mesh,
	const fileName&,
	const triSurface&,
	const triSurfaceSearch& querySurf,
	const pointField& outsidePts,

	const label nCutLayers,

	labelHashSet& inside,
	labelHashSet& outside,
	labelHashSet& cut
)
{
	// Construct search engine on mesh
	meshSearch queryMesh(mesh, true);


	// Check all 'outside' points
	forAll(outsidePts, outsideI)
	{
		const point& outsidePoint = outsidePts[outsideI];

		// Find cell point is in. Linear search.
		if (queryMesh.findCell(outsidePoint, -1, false) == -1)
		{
			FatalErrorIn
			(
				"surfaceSets::getSurfaceSets"
				"(const polyMesh&, const fileName&, const triSurface&"
				", const triSurfaceSearch&, const pointField&"
				", labelHashSet&, labelHashSet&, labelHashSet&)"
			)   << "outsidePoint " << outsidePoint
				<< " is not inside any cell"
				<< exit(FatalError);
		}
	}

	// Cut faces with surface and classify cells
	cellClassification cellType
	(
		mesh,
		queryMesh,
		querySurf,
		outsidePts
	);

	if (nCutLayers > 0)
	{
		// Trim cutCells so they are max nCutLayers away (as seen in point-cell
		// walk) from outside cells.
		cellType.trimCutCells
		(
			nCutLayers,
			cellClassification::OUTSIDE,
			cellClassification::INSIDE
		);
	}

	forAll(cellType, cellI)
	{
		label cType = cellType[cellI];

		if (cType == cellClassification::CUT)
		{
			cut.insert(cellI);
		}
		else if (cType == cellClassification::INSIDE)
		{
			inside.insert(cellI);
		}
		else if (cType == cellClassification::OUTSIDE)
		{
			outside.insert(cellI);
		}
	}
}


Foam::labelHashSet Foam::surfaceSets::getHangingCells
(
	const primitiveMesh& mesh,
	const labelHashSet& internalCells
)
{
	const cellList& cells = mesh.cells();
	const faceList& faces = mesh.faces();


	// Divide points into
	// -referenced by internal only
	// -referenced by outside only
	// -mixed

	List<pointStatus> pointSide(mesh.nPoints(), NOTSET);

	for (label cellI = 0; cellI < mesh.nCells(); cellI++)
	{
		if (internalCells.found(cellI))
		{
			// Inside cell. Mark all vertices seen from this cell.
			const labelList& cFaces = cells[cellI];

			forAll(cFaces, cFaceI)
			{
				const face& f = faces[cFaces[cFaceI]];

				forAll(f, fp)
				{
					label pointI = f[fp];

					if (pointSide[pointI] == NOTSET)
					{
					    pointSide[pointI] = INSIDE;
					}
					else if (pointSide[pointI] == OUTSIDE)
					{
					    pointSide[pointI] = MIXED;
					}
					else
					{
					    // mixed or inside stay same
					}
				}
			}
		}
		else
		{
			// Outside cell
			const labelList& cFaces = cells[cellI];

			forAll(cFaces, cFaceI)
			{
				const face& f = faces[cFaces[cFaceI]];

				forAll(f, fp)
				{
					label pointI = f[fp];

					if (pointSide[pointI] == NOTSET)
					{
					    pointSide[pointI] = OUTSIDE;
					}
					else if (pointSide[pointI] == INSIDE)
					{
					    pointSide[pointI] = MIXED;
					}
					else
					{
					    // mixed or outside stay same
					}
				}
			}
		}
	}


	//OFstream mixedStr("mixed.obj");
	//
	//forAll(pointSide, pointI)
	//{
	//    if (pointSide[pointI] == MIXED)
	//    {
	//        const point& pt = points[pointI];
	//
	//        mixedStr << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z()
	//            << endl;
	//    }
	//}


	// Determine cells using mixed points only

	labelHashSet mixedOnlyCells(internalCells.size());

	for
	(
		labelHashSet::const_iterator iter = internalCells.begin();
		iter != internalCells.end();
		++iter
	)
	{
		label cellI = iter.key();

		const cell& cFaces = cells[cellI];

		label usesMixedOnly = true;

		forAll(cFaces, i)
		{
			const face& f = faces[cFaces[i]];

			forAll(f, fp)
			{
				if (pointSide[f[fp]] != MIXED)
				{
					usesMixedOnly = false;
					break;
				}
			}

			if (!usesMixedOnly)
			{
				break;
			}
		}
		if (usesMixedOnly)
		{
			mixedOnlyCells.insert(cellI);
		}
	}

	return mixedOnlyCells;
}


//void Foam::surfaceSets::writeSurfaceSets
//(
//    const polyMesh& mesh,
//    const fileName& surfName,
//    const triSurface& surf,
//    const triSurfaceSearch& querySurf,
//    const pointField& outsidePts,
//    const scalar edgeFactor
//)
//{
//    // Cellsets for inside/outside determination
//    cellSet rawInside(mesh, "rawInside", mesh.nCells()/10);
//    cellSet rawOutside(mesh, "rawOutside", mesh.nCells()/10);
//    cellSet cutCells(mesh, "cutCells", mesh.nCells()/10);
//
//    // Get inside/outside/cut cells
//    getSurfaceSets
//    (
//        mesh,
//        surfName,
//        surf,
//        querySurf,
//        outsidePts,
//
//        rawInside,
//        rawOutside,
//        cutCells
//    );
//
//
//    Pout<< "rawInside:" << rawInside.size() << endl;
//
//    label nRemoved;
//    do
//    {
//        nRemoved = removeHangingCells(mesh, querySurf, rawInside);
//
//        Pout<< nl
//            << "Removed " << nRemoved
//            << " rawInside cells that have all their points on the outside"
//            << endl;
//    }
//    while (nRemoved != 0);
//
//    Pout<< "Writing inside cells (" << rawInside.size() << ") to cellSet "
//        << rawInside.instance()/rawInside.local()/rawInside.name()
//        << endl << endl;
//    rawInside.write();
//
//
//
//    // Select outside cells
//    surfaceToCell outsideSource
//    (
//        mesh,
//        surfName,
//        surf,
//        querySurf,
//        outsidePts,
//        false,          // includeCut
//        false,          // includeInside
//        true,           // includeOutside
//        -GREAT,         // nearDist
//        -GREAT          // curvature
//    );
//
//    outsideSource.applyToSet(topoSetSource::NEW, rawOutside);
//
//    Pout<< "rawOutside:" << rawOutside.size() << endl;
//
//    do
//    {
//        nRemoved = removeHangingCells(mesh, querySurf, rawOutside);
//
//        Pout<< nl
//            << "Removed " << nRemoved
//            << " rawOutside cells that have all their points on the outside"
//            << endl;
//    }
//    while (nRemoved != 0);
//
//    Pout<< "Writing outside cells (" << rawOutside.size() << ") to cellSet "
//        << rawOutside.instance()/rawOutside.local()/rawOutside.name()
//        << endl << endl;
//    rawOutside.write();
//
//
//    // Select cut cells by negating inside and outside set.
//    cutCells.invert(mesh.nCells());
//
//    cellToCell deleteInsideSource(mesh, rawInside.name());
//
//    deleteInsideSource.applyToSet(topoSetSource::DELETE, cutCells);
//    Pout<< "Writing cut cells (" << cutCells.size() << ") to cellSet "
//        << cutCells.instance()/cutCells.local()/cutCells.name()
//        << endl << endl;
//    cutCells.write();
//
//
//    //
//    // Remove cells with points too close to surface.
//    //
//
//
//    // Get all points in cutCells.
//    pointSet cutPoints(mesh, "cutPoints", 4*cutCells.size());
//    cellToPoint cutSource(mesh, "cutCells", cellToPoint::ALL);
//    cutSource.applyToSet(topoSetSource::NEW, cutPoints);
//
//    // Get all points that are too close to surface.
//    pointSet nearPoints(mesh, "nearPoints", cutPoints.size());
//
//    getNearPoints
//    (
//        mesh,
//        surf,
//        querySurf,
//        edgeFactor,
//        cutPoints,
//        nearPoints
//    );
//
//    Pout<< nl
//        << "Selected " << nearPoints.size()
//        << " points that are closer than " << edgeFactor
//        << " times the local minimum lengthscale to the surface"
//        << nl << endl;
//
//
//    // Remove cells that use any of the points in nearPoints
//    // from the inside and outsideCells.
//    nearPoints.write();
//    pointToCell pToCell(mesh, nearPoints.name(), pointToCell::ANY);
//
//
//
//    // Start off from copy of rawInside, rawOutside
//    cellSet inside(mesh, "inside", rawInside);
//    cellSet outside(mesh, "outside", rawOutside);
//
//    pToCell.applyToSet(topoSetSource::DELETE, inside);
//    pToCell.applyToSet(topoSetSource::DELETE, outside);
//
//    Pout<< nl
//        << "Removed " << rawInside.size() - inside.size()
//        << " inside cells that are too close to the surface" << endl;
//
//    Pout<< nl
//        << "Removed " << rawOutside.size() - outside.size()
//        << " inside cells that are too close to the surface" << nl << endl;
//
//
//
//    //
//    // Remove cells with one or no faces on rest of cellSet. Note: Problem is
//    // not these cells an sich but rather that all of their vertices will be
//    // outside vertices and thus projected onto surface. Which might (if they
//    // project onto same surface) result in flattened cells.
//    //
//
//    do
//    {
//        nRemoved = removeHangingCells(mesh, querySurf, inside);
//
//        Pout<< nl
//            << "Removed " << nRemoved
//            << " inside cells that have all their points on the outside"
//            << endl;
//    }
//    while (nRemoved != 0);
//    do
//    {
//        nRemoved = removeHangingCells(mesh, querySurf, outside);
//
//        Pout<< nl
//            << "Removed " << nRemoved
//            << " outside cells that have all their points on the inside"
//            << endl;
//    }
//    while (nRemoved != 0);
//
//
//    //
//    // Write
//    //
//
//
//    Pout<< "Writing inside cells (" << inside.size() << ") to cellSet "
//        << inside.instance()/inside.local()/inside.name()
//        << endl << endl;
//
//    inside.write();
//
//    Pout<< "Writing outside cells (" << outside.size() << ") to cellSet "
//        << outside.instance()/outside.local()/outside.name()
//        << endl << endl;
//
//    outside.write();
//}


// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //


// * * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * //


// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //


// ************************************************************************* //
